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A cut-cell approach to Computational Fluid Dynamics (CFD) that utilizes the median 
dual of a tetrahedral background grid is described. The discrete adjoint is also calculated, 
which permits adaptation based on improving the calculation of a specified output (off-body 
pressure signature) in supersonic inviscid flow. These predicted signatures are compared 
to wind tunnel measurements on and off the configuration centerline 10 body lengths below 
the model to validate the method for sonic boom prediction. Accurate mid-field sonic boom 
pressure signatures are calculated with the Euler equations without the use of hybrid grid 
or signature propagation methods. Highly-refined, shock-aligned anisotropic grids were 
produced by this method from coarse isotropic grids created without prior knowledge of 
shock locations. A heuristic reconstruction limiter provided stable flow and adjoint solu- 
tion schemes while producing similar signatures to Barth-Jespersen and Venkatakrishnan 
limiters. The use of cut-cells with an output-based adaptive scheme completely automated 
this accurate prediction capability after a triangular mesh is generated for the cut surface. 

This automation drastically reduces the manual intervention required by existing methods. 

I. Introduction 

The acceptance of an aircraft’s sonic boom to the general population 
is a requirement for supersonic flights over land and therefore the com- 
mercial viability of a supersonic transport. Predicting how sonic boom 
signatures are perceived is a challenging task that requires the prediction 
of the signature on the ground. This is a task complicated by the long 
propagation distances, atmosphere variations, and the Earth’s turbulent 
boundary layer. A detailed review of the history and state-of-the-art of 
sonic boom modeling is provided by Plotkin. 1 

The propagation of a sonic boom is often separated into three logical 
stages or regions, depicted in Fig. 1, to facilitate analysis. 2 The near- 
field is a region near the aircraft, where shocks are formed and strongly 
influenced by nonlinear phenomena such as shock-shock interaction, shock 
curvature, and cross flow. Higher pressure portions of the signature travel 
faster than lower pressure portions of the signature because of variations 
in the local speed of sound. This slight speed difference causes the shocks 
to deform by elongating and coalescing in the mid-field. The signature 
is also refracted by variations in the atmospheric speed of sound. The 
signature has become primarily an N-wave in the far-held as a result of 
this distortion. The boundaries of these regions are case specific. 
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Whitham 3,4 provides analytic solutions for the signal distortion of slender axisymmetric projectiles. 
Boom propagation has been implemented in a number of computer programs. 5,6 Unfortunately, these boom 
propagation methods are not directly applicable to complex aircraft geometries. Page and Plotkin 7 applied 
the multipoles of George 8 to combine CFD near-held calculation with mid- and far-held boom propagation. 
This CFD matching multipole propagation technique has been revisited by Rallabhandi and Mavris. 9 

CFD codes have difficulty propagating the relatively weak pressure signatures of a sonic boom to distances 
beyond the near-held region, where these boom propagation methods are valid. This problem is more acute 
for unstructured grid methods that are often employed to capture the geometrical complexity of the model, 
especially if the grids are not aligned with the shocks. To improve alignment, isotropic unstructured grids are 
stretched to align the tetrahedra with the free stream Mach angle to improve signal propagation for initial 
grids. 10 This alignment issue has also given rise to hybrid methods 11-13 where near-body unstructured 
grid solutions are interpolated to shock-aligned structured grid methods to increase accuracy. The hybrid 
methods are hindered by the interpolation process, so adaptive grid methods 14, 15 are employed to improve 
the accuracy of unstructured grid methods for long propagation distances. These previous adaptive methods 
have used only primal solution information (Mach and density) to drive the adaptive process. 

Adaptive grid methods are designed to specify a resolution request that equidistributes and minimizes 
error estimates. The improved resolution request is commonly based on local error estimates. 16-18 Uniformly 
reducing the errors associated with all local-error sources of the flow may not be optimal from an engineering 
context, where calculating an output functional (i.e. , boom signature) may be of greater concern. An 
alternative method is to estimate the error in the calculation of a specified engineering output functional. 19-22 
Output error indicators utilize the dual or adjoint solution of an output functional to account for the impact 
of local error as well as the transport of these local errors throughout the problem domain to improve 
the calculation of that output functional. This output-adaptive approach has been applied to sonic boom 
prediction in 2D with discontinuous Galerkin 23, 24 and Cartesian methods. 22 

In this work, anisotropic output-based adaptation to improve an off-body pressure integral is applied 
to 3D sonic boom prediction. Anisotropically adapted tetrahedral background grids with cut-cells provide 
extremely robust adaptation mechanics, enabling the automated application of anisotropic output-based 
adaptation to non-trivial 3D sonic boom problems for the first time. This allows for the entire signature to 
be calculated or the pressure integral can be restricted to a specific region of interest. The adjoint solution 
can also provide engineering intuition with a rigorous foundation for design sensitivities 25 and discretization 
error estimates. 21 

Cut-cell methods with Cartesian background grids 26-29 have been very successful for Euler simulations. 
The regular structure of the Cartesian background grid permits extremely efficient solution schemes. Carte- 
sian background grids have the capability to only provide anisotropic resolution in the Cartesian directions. 29 
Simplex meshes have the ability to stretch the triangular and tetrahedral elements in arbitrary directions. 
This permits the efficient representation of anisotropic features (i.e., shocks). The cut-cell method is also 
applicable to simplex meshes. 30-32 When the constraint of providing a body-fitted grid is removed, the grid 
adaptation task becomes much simpler. The complexities of adaptation on curved domain boundaries 33 is 
eliminated and robustness is dramatically increased. 

The 3D anisotropic output-adaptive method for sonic boom prediction has previously been applied to 
body-fitted grids. 34,35 This work addresses two of the weaknesses uncovered in the previous work: adjoint 
iterative convergence for reconstruction limiters and the robustness of body- fitted adaptive mechanics. The 
current method is validated by comparison to wind tunnel data for representative configurations. 

II. Cut-Cell Determination 

To introduce the 3D cut-cell method a simple 2D example is presented. The 
control volumes used by the flow solver are the median duals of a triangular mesh. 

These median duals are constructed by gathering the three dual faces that are inside 
each primal triangle. Figure 2 shows the three dual faces (dashed lines) for a triangle, 
which each connect the triangle center to one of the triangle side midpoints. A 
primal triangular grid is shown in Fig. 3(a). The median duals of this triangular 
grid are shown in Fig. 3(b). The airfoil geometry is a diamond airfoil, shown with 
the uncut median dual background grid in Fig. 3(c). The airfoil geometry is Boolean 
subtracted from this background grid removing the portion of the background grid 



Figure 2. Three dual 
faces (dashed lines) as- 
sociated with a single 
triangle (solid lines). 
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that is external to the flow domain, Fig. 3(d). 

In the 3D case, the domain of the simulation is constructed by Boolean subtraction of a manifold tri- 
angular boundary representation from a background grid. This triangulation can come from many sources. 
Figure 4(b) is a triangular surface grid of a cylinder constructed on a CAD solid. 36 



(c) Median dual grid with geometry. 


(d) Resulting cut-cell grid. 


Figure 3. Cut cell illustration of an diamond airfoil in 2D. 



(a) Median dual. 


(b) Cylindrical cutting surface. 


(c) Resulting Polyhedra. 


Figure 4. Cylindrical cutting surface and median dual. 


The background grid contains closed simplicial polytope control vol- 
umes. In 3D, these polyhedra are the median duals of a tetrahedral grid. 
The 3D median dual about a single primal node is shown in Fig. 4(a). 
Just as in the 2D case, this dual control volume may not be convex. Each 
dual polyhedra of a tetrahedral grid contains 0(100) triangles. Figure 
5 illustrates the two dual triangles associated with an edge of a primal 
tetrahedron. There are six edges in a tetrahedron, which contains a total 
of 12 triangular dual faces shared by the duals at each of its 4 nodes. For 
robustness and a decrease in execution time and memory usage, a trian- 
gular dual face is only represented once in the intersection procedure and 
shared by the two adjacent control volumes. 



Figure 5. Two dual triangles associ- 
ated with a single tetrahedral edge. 
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Aftosmis, Berger, and Melton 28 demonstrate a method to perform the 
Boolean subtraction of two manifold triangular polyhedra (surface grid and each of the background grid 
control volumes) with a series of triangle-triangle intersections. Their methodology was modified to use only 
floating point arithmetic by Park. 3 ' The result of this subtraction is shown in Fig. 4(c). 


III. Flow and Adjoint Solvers 


Fully Unstructured Navier-Stokes Three-Dimensional (FUN3D) is a suite of codes for finite-volume 
CFD . 38 The FUN3D website, http : //f un3d . larc . nasa . gov, contains the user manual and an extensive list 
of references. FUN3D is able to solve incompressible, Euler, and Reynolds-averaged Navier-Stokes (RANS) 
flow equations, either tightly or loosely coupled to a turbulence model. The Euler equations are used in this 
study. Domain decomposition is employed to fully exploit the distributed memory of a cluster of computers 
to increase problem size and reduce the execution time of the simulation process. 

The Euler equations are 
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where p is density, u, v , and iu are velocity, E is total energy per unit volume, and p is pressure. These 
quantities are related by the ideal gas relation, 
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with the specific heat ratio 7 = 1.4 for air. 

The divergence theorem is applied over a set of control volumes to produce a finite-volume scheme, 
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where T; are the boundaries of the control volumes with volume V) and n is an outward pointing normal. 
The average of Q in each control volume is Q*. The flux integration is approximated as, 
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where Ri is the discrete residual for control volume i, the summation is over the faces of the control volume . 38 
The van Leer 39 approximate Riemann solver Ti is utilized to compute the flux from the primitive states, 
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Q = 
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at the borders of the neighboring control volumes, q r f and qi j. These face values are reconstructed from cell 
averages (the reconstruction method is described below). The discrete equations are established simultane- 
ously for each control volume, 

F^ + R(Q) = °, (7) 

which makes the discrete solution vector Q £ 1Z 5N , discrete residual vector R £ 1Z 5N , and V = diag(F)), 
where N is the number of control volumes. The flux integration scheme (including face state reconstruction 
from cell averages) is detailed in the following sections. 

A backward Euler solution update scheme is employed with a variable pseudo-time step . 38 An approxi- 
mate nearest neighbor linearization is utilized to reduce the memory required for the implicit point-iterative 
method. After the flow solution is known, the discrete adjoint equations 25,40 are solved to complete the 
dual problem. The linear adjoint equations are solved with a dual-consistent time-marching method . 41,42 
The dual-consistent solution method guarantees that the adjoint equations will have the same asymptotic 
convergence rate as the flow equations. 
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III. A. Inviscid Flux Integration 


The existing FUN 3D body- fitted approach lumps the median dual pieces to a single effective area and normal 
direction for each edge they surround. 43 After lumping, all of the inviscid terms are calculated with a loop 
over edges, which is computationally efficient. Conserved states Q , used in the time advancement scheme, 
are converted to primitive states q for face state reconstruction. The primitive state is extrapolated from 
the nodes to establish the primitive state at these lumped faces Qf using the gradients Vq = [q x ,q y ,q z ], face 
center Xf, and node xq, 

qf = Qo + Vq (x f - Zq), (8) 


for the unlimited scheme. The gradient is reconstructed from the cell averaged state qo with a weighted least 
squares method, see Park 37 for details. For the case of supersonic flow, a limiting function is used to reduce 
the gradient contribution to the reconstruction (see Subsection III.B). 




(c) Cut-cell boundary integra- 
tion. 


Figure 6. Dual control volumes, in 2D. 


At the completion of cut-cell preprocessing, the dual polyhedra can be classified into three groups: uncut 
active duals interior to the computational domain, cut duals, and inactive uncut duals exterior to the 
computational domain. The state is stored at each node in the primal grid, Fig. 6(a) filled circle. All nodes 
that correspond to dual polyhedra that have been cut or are inactive are removed. A new degree of freedom 
is inserted at each cut dual polyhedra centroid, Fig. 6(b) filled circle. Multiple degrees of freedom are added 
when a polyhedra is split into multiple distinct regions by the cut surface. 

Once a dual control volume is cut, the approximation that the state is centered at the primal node is 
removed and the state becomes centered at the control volume centroid. This results in a discontinuous 
change in location once the control volume is infinitesimally cut. This discontinuous behavior may cause 
difficulties for shape sensitivities and design. Removing this issue remains a topic for future work, but may 
be addressed by computing the uncut cell centroids. 

The median dual triangles that surround any edges that involve a cut cell are explicitly represented and 
employed in flux integration. Edges that involve uncut active duals utilize the lumped effective areas and 
normals of the body-fitted scheme. Cut-cell flux integration requires more work than the body-fitted scheme 
because there are multiple triangles separating the two control volumes that would be approximated as a 
single flux evaluation in the body-fitted scheme. It also requires more memory to store the extra triangles 
that would be approximated as a single effective area. The cut cells are a minority of the control volumes 
for a typical case, so the additional expense of utilizing cut cells does not dominate the execution time or 
storage. 

The body-fitted, node-based scheme stores the state on the boundary of the domain. The state is 
interpolated between adjacent boundary nodes to integrate the boundary flux. The boundary flux is described 
in Section IV. For cut-cell boundary flux integration, the state is extrapolated with the reconstructed 
gradients from the cell centroids to the boundary face, 

qbf = qo + Vq (x bf - £ 0 ). (9) 

The boundary of the cut cell from Fig. 6(b) is shown in Fig. 6(c) to illustrate the reconstruction of the 
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boundary state qb f for cut-cell boundary flux integration. The qb f is extrapolated to each cut surface piece 
for integration. 


III.B. Reconstruction Limiting 

Barth and Jespersen 43 introduced limits on an unstructured grid reconstruction scheme to maintain mono- 
tonicity. Face reconstruction using a limiter of this form is 

Qf =% + $ Vq(xf -x 0 ), (10) 


where the diagonal matrix limiting function is computed in each control volume. The same $ is employed 
in all face reconstructions for a given control volume. This type of limiter can compromise the convergence of 
the flow solution and therefore a dual-consistent adjoint solver. 44,45 Venkatakrishnan 46 studied this limiter 
in its original form as well as with the limiter function held constant (i.e. frozen) after iterative convergence 
stalls. He proposed a new limiter to improve convergence, but both the frozen scheme and new limiter can 
result in stalled convergence. The Venkatakrishnan limiter is not monotone, it permits under- and over- 
shoots. Frozen limiters are approximations that impede error estimation, output-based adaptation, adjoint 
iterative convergence, and design sensitivities. 22,34,4 ' 

In this study, the limiter will be used in the context of an output adaptive scheme that requires the 
adjoint solution. An exact linearization and steady iterative convergence of the flow and adjoint solvers 
is paramount to the robustness of the adaptive scheme. This iterative convergence is so critical that the 
accuracy of the limited scheme will be sacrificed; accuracy will be regained with adaptive grid refinement 
and alignment. A heuristic edge-based limiter 48 is utilized to improve the convergence of the flow solver 
while providing the exact linearization required for adjoint convergence. Concessions are made to improve 
iterative convergence; it is not total variation diminishing (TVD) or linearity preserving. 

The heuristic limiter was developed 48 by examining its effect on shock capturing for regular and irregular 
grids and empirically adjusting its formulation to increase the width of shocks. It is a scalar limiting 
function <j> that considers only the cell-averaged values of pressure and their reconstructed gradients in the 
cells adjacent to the face being reconstructed. Face reconstruction using a limiter of this form is 

q f = q + (/)Vq(x f - x 0 ), (11) 

where the scalar limiting function <f> is computed for each face /. The same is used for the left qif and 
right q r j face reconstructions. 



Figure 7. Edge and face geometry. 


The basic concept employed in this heuristic limiter is to reduce the reconstruction gradient in locations 
where the pressure gradients are large relative to pressure. This clearly could result in limiting in regions for 
which the solution varies linearly (though with large magnitude), however, in combination with adaptation 
the proposed limiter has been found robust and accurate. The specific form of the limiter relies on a measure 
of the change in the pressure, Sp. To form Sp , the reconstructed gradient of pressure for the control volumes 
on the right and left of the face (Fig. 7), Vpi and Vp 2 , are used with the right and left extrapolation vectors 
to the face, Sx\ and Sx 2 , 


6p = 


Sx iyPK - Sx 2 X ^P2 X 
8x\ Vpi — Sx 2 V»2 
Sxi z Vpi z - Sx 2z Vp 2z 


(12) 


This sensor is active for linear functions and does not specifically penalize extrema. The gradient recon- 
struction is reduced where the the Sp sensor is large with the intention of spreading the detected jump over 
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a number of control volumes. Adaptation will be employed to narrow the width of the discontinuity. The 
tanh function is employed to smooth the combined nondimensional pressure jump ratio, 


^heuristic — 1 tanh 


Sp 

min(pi,p 2 ) 


(13) 


and restricts the limiter to the range (0,1]. A tanh function is employed to provide a smoothly varying 
and differentiable function that enables residual convergence that can be impeded by a non-smooth limiting 
function. This limiter is active (to some degree) in all regions with pressure variations, so it will not switch 
on and off intermittently during iterative convergence. The design accuracy of the limited scheme is therefore 
below second-order. The limiter is more active when the pressure variation is significant as compared to the 
local pressure. 

The cut cells require pressure extrapolated to the boundaries to compute boundary fluxes. This recon- 
struction requires limiting to prevent unrealizable face states and must be smoothly differentiable to facilitate 
iterative convergence, 


Spd = 


Sx x V Px 

SXyVPy , 

5x z Vp z 


(14) 


^extrapolation 1 tanh 


(15) 


The extrapolation limiter is formulated to mimic the interior face limiter using only the data from the cell 
adjacent to the boundary. 


IV. Boundary Conditions 

The boundary conditions are imposed weakly through the fluxes. The tangential flow boundary condition 
is implemented with zero velocity normal to the boundary, resulting in the flux, 


F t 


tangential 


o 

pn x 

pn y , 
pn z 
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(16) 


where p is interpolated along or extrapolated to the boundary. The supersonic outflow boundary condition 
uses the interior state to form the boundary flux. The supersonic inflow boundary condition uses the free 
stream state poo, Uqq, i>oo, u>oo, and p^ to form the boundary flux. 

The inviscid flow model breaks down at a sharp corner where separation would occur in a physical flow. 
In the supersonic flow simulations performed in this work, this problem presents itself at blunt trailing 
edges. To avoid this problem in these regions, a transpiration boundary condition is specified manually. 
This boundary condition applies free stream velocity state, ttoo, Uqo, and with a density and pressure of 
p = 0.3poo and p = 0.3poo. This level of density and pressure is empirically established by examining the 
solution of a backward facing step with the tangential boundary condition. 


V. Output-Based Adaptation 

Venditti 21 describes an output-based error estimation and adaptation scheme. To formulate the error 
estimate, an embedded grid is required. Constructing the entire embedded grid can be infeasible for large 
3D grids and has prevented the use of adjoint error estimation techniques for industrial-sized problems even 
with a parallel implementation. 34 While the embedded grid can be formed in sections, this increases the 
error estimation scheme complexity. Forming a portion or the entire embedded grid is also complicated by 
the need to respect curved boundaries and recompute the intersection tests of cut cells. These difficulties 
have motivated the desire to employ only the current grid in the error estimation procedure. A procedure 
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is described that obtains an indicator for output adaptation with the current grid, but does not provide a 
functional error correction. 

Park 3 ' provides a derivation of the single-grid adaptive indicator by placing it in the context of the 
embedded grid approach of Venditti. In the interest of brevity, the single-grid error estimation and adaptive 
indicator is provided without a derivation, 


[-^single ] k = 2 5Z { 


[A - A] i)(e [i?(Q)]i iK |} 


(17) 


It has the same pieces as the Venditti error indicator, where the five conservation equations are contracted 
by the summation over i. The vector / s i ng i e £ 1Z N has a single value for each grid control volume k. The 
A and Q higher-order reconstructions and the A, and Q lower-order reconstructions on the current grid are 
described in Park. 37 The original residual operators are utilized and A and Q are constructed to make i? A ( A) 
and R{Q) reliable adaptive indicators. The A and Q reconstructions are formed with a fit of quadratic 
functions to cell averaged states and their gradients. The difference between the () and () reconstructions is 
intended to provide adequate guidance for the relative distribution of error, not a sharp bound on error. 

Venditti 21 provides a procedure to calculate a new grid spacing request h from the adaptive indicator 
I K and an error tolerance tolo. The adjoint adaptation parameter was also incorporated into an anisotropic 
Hessian-based framework. This combined approach sets the anisotropy of mesh elements by using the Mach 
Hessian, and it scales the element size so that the tightest spacing is dictated by the adjoint adaptation 
parameter. The metric-based grid mechanics is described by Park and Darmofal. 37,49 It is fully parallelized 
and interfaces directly with the error estimation process. 


VI. Application to Sonic Boom Prediction 


The parallel metric-based adaptation algorithm is applied to sonic boom prediction. More detailed 
information on these cases, including timing information, is available in Park. 3 ' The output for the adaptive 
procedure of the integral of quadratic pressure deviation over a surface s in the domain, 


/ = 


1 

A s 


f p- Poo 
V Poo 


2 

els, 


(18) 


where A s is the area of the integration surface. This focuses the adaptation on improving the calculation 
of pressure near this surface. Previous applications have been performed with the integral of pressure 
deviation. 34,35 However, the square of this deviation has been shown to produce more accurate signatures 
with less control volumes. 22 A cylindrical integration surface, aligned to the a;-axis, is employed. The extent 
of the integration surface can be optionally restricted to the interior of a box to focus on a subsection of the 
cylinder. 

The current output-based adaptation approach is validated with wind tunnel measurements. Wind 
tunnel testing of sonic boom configurations is a challenging task. Wind tunnel models are typically small 
to obtain pressure signatures a relatively large distance from the model within the finite size of tunnel test 
sections. Carlson and Morris 50 present some of difficulties inherent in wind tunnel testing of these small 
models including extraneous variations in pressure larger than the signals measured. A test apparatus and 
procedure that mitigates the extraneous spatial and temporal distortions is described. Morgenstern 51 also 
documents variations in ambient static pressure wind tunnel measurements that are of the same magnitude 
as the desired signature measurement. 

VI. A. Cone-Cylinder Configuration 

A double cone geometry, denoted “Model 8” in a 1965 wind tunnel report, 52 
is shown in Fig. 8 with a shaded triangular surface grid. This same case was 
employed to evaluate 34 and then validate 35 a parallel adaptive body-fitted grid 
approach. This configuration has also been used by other researchers to eval- 
uate their signature prediction techniques. 12, 15 The pressure integral output 
function was defined as a cylinder, six body lengths in radius, centered about 
the geometry axis. The cylinder is clipped forward of 3 body lengths behind the 



Figure 8. Cone-Cylinder Ge- 
ometry. 
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nose, aft of 9 body lengths behind the nose, and outside of 0.1 body lengths off the centerline to focus only 
on the region where wind tunnel data is available. The surface grid is Boolean subtracted from a 9-degree 
wedge-shaped background tetrahedral volume grid. A symmetry plane of the volume grid is shown in Figure 
9. The initial grid (4,000 control volumes) was created with no prior knowledge of where the shocks would 
propagate through the domain, Fig. 9(a). The free stream Mach number is 1.26 and the heuristic limiter is 
employed during adaptation. The parallel execution scheme used 32 partitions, and the 17 th adapted grid 
(7,500,000 control volumes) is shown in Fig. 9(b). The shocks have been implicitly targeted and refined to 
propagate the signal to the pressure integral surface. The anisotropy of the grid, based on the Mach Hessian, 
is clearly evident. This anisotropy reduces the number of required control volumes. 



(a) Initial grid. 


(b) Output adapted grid. 


Figure 9. Symmetry planes of initial and output adapted double cone 3D volume grids. 


The adaptation history of the pressure integral and its remaining error estimation is shown in Fig. 10. 
The remaining error estimate is given by Eq. (17). Error is underestimated on the initial few adapted grids 
before the shocks are propagated to the integration surface. Once this connection is established, the error is 
reduced. 

The adaptation history of the pressure signature extracted at 6 body lengths is shown in Fig. 11. The 
circular symbols are digitized from a wind tunnel report. 52 The solid line is the final, adapted signature. 
The signal is absent on the original coarse grid. The extrema of the pressure signal start to form and grow. 
The inflection points at x/l = 2.3 are the last part of the signal to form. The over- and under-shoots of the 
signal intensify on the final few adapted grids as the grid-shock alignment improves. 

The final adapted grid is simulated with the Venkatakrishnan, heuristic, and Barth- Jespersen limiters 
in Fig. 12. The Venkatakrishnan limiter has similar over- and under-shoots to the heuristic limiter. The 
Barth-Jespersen limiter produces a signature without over- and under-shoots. All of these limiters have very 
similar signatures except at the discontinuities. 

VI. B. Delta-Wing Configuration 

A delta wing body, denoted “Model 4” in a 1973 wind tunnel report, 53 is shown in Fig. 13 with the triangular 
surface grid. This case has also been used by other researchers to evaluate their techniques. 54-58 The delta 
wing body was simulated at Mach 1.68 and 0.0 deg angle of attack. The cylindrical integration surface 
has a radius of 3.6 body lengths and is clipped fore ( x/l = 4.3) and aft (x/l = 6.6) to restrict its extent 
to the location of the wind tunnel data. These clipping locations are 0.6 body lengths ahead of and 1.7 
body lengths behind of the free stream Mach cone emanating from the nose. The width of the integration 
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(a) Pressure integral. 


(b) Pressure integral remaining error. 


Figure 10. Model 8 pressure integral and uncertainty convergence at 6 body lengths. 



Figure 11. Model 8 pressure signature adaptation history at 6 body lengths. 



Figure 12. 


Model 8 final adapted pressure signature at 6 body lengths for various limiters. 
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Figure 13. Delta Wing Geometry. 


surface is clipped to half a wingspan. Figure 14 shows the initial (2,800 control volumes) and final adapted 
(4,900,000 control volumes) grids colored with pressure. The background grid exit plane and symmetry 
plane intersected with the model surface grid is shown. The initial grid is extremely coarse. The model 
is contained in a small number of initial background grid control volumes. The final adapted grid is well 
aligned with the propagated signal. The refinement region coarsens very rapidly away from the symmetry 
plane (the integration surface is only a half wingspan wide). 



(a) initial grid. 


(b) output adapted grid. 


Figure 14. Exit and symmetry planes of initial and output adapted delta wing body background grids. 


The adaptation history of the pressure integral and remaining error estimation is shown in Fig. 15. Error 
is initially under predicted for this extremely coarse initial grid. The robustness of the adaptive procedure 
is illustrated by the use of this coarse initial grid. The remaining error estimate improves as the shocks 
are resolved and propagated to the integration surface on grids smaller than 100,000 control volumes. This 
remaining error estimate steadily decreases on grids larger than 100,000 control volumes as the signature is 
refined. The grid has grown in size by three orders of magnitude. 

Figure 16 shows the pressure signatures for the series of grids employed during adaptation. The pressure 
signature is not apparent on the initial grid. The signals grow as a result of adaptation and the over- and 
under-shoots increase on the final few grids. The heuristic limiter is employed during adaptation. Figure 
17 shows the final adapted grid pressure signatures for the Venkatakrishnan, heuristic, and Barth-Jespersen 
limiters. The signatures of all three limiters are very similar, except near discontinuities. The Barth- 
Jespersen limiter eliminates the over- and under-shoots of the bow and tail shocks. The difference between 
the Venkatakrishnan and heuristic limiter signatures is very small. 
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(a) Pressure integral. (b) Pressure integral remaining error. 

Figure 15. Delta wing body pressure integral and uncertainty convergence at 3.6 body lengths. 



Figure 16. Delta wing body pressure signature adaptation history at 3.6 body lengths. 



Figure 17. Delta wing body final adapted pressure signature at 3.6 body lengths for various limiters. 
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Figure 18. Low-boom wing body geometry. 


VI. C. Low-Boom Configuration 

The Straight-Line Segmented Leading Edge (SLSLE) low-boom configuration, Fig. 18, is described by Mack 
and Kuhn. 59,60 These reports provide wind tunnel data from two tests, performed at the Langley Research 
Center Unitary Plan Wind Tunnel Facility 59 and the John Glenn Research Center 10 x 10 ft Wind Tunnel 
Facility. 60 The test condition is Mach 2.0. The model surface geometry for the aircraft includes the model 
and sting incidence that provides the wind tunnel lift coefficient Cl = 0.08309 so the angle of attack is 
zero. 11 The configuration has a finite thickness trailing edge, which was modeled with a transpiration 
boundary condition to prevent a strong inviscid supersonic corner flow expansion. Preliminary body- fitted 
results for this configuration 34 extended the blunt trailing edge to a sharp trailing edge to avoid the strong 
supersonic expansion. 

The original symmetry plane and cut surface grid colored with pressure is shown in Fig. 19(a). A linear 
distribution of pressure is shown in each control volume, resulting in a discontinuous pressure distribution on 
the surface. The initial background grid is isotropic. The final adapted symmetry plane and cut surface grid 
colored with pressure is shown in Fig. 19(b). The initial background grid contains 40,000 control volumes, 
and the final adapted background grid contains 5,700,000 control volumes. The anisotropy of the adapted 
grid has been established with the Mach Hessian, aligning the grid with the shocks. 



(a) Initial grid. (b) Final grid. 

Figure 19. SLSLE surface grid colored with pressure. 


The initial, Fig. 20(a), and final, Fig. 20(b), grid integration surfaces are colored with pressure deviation 
from free stream. These cylindrical integration surfaces are used to compute the output, which drives the 
adaptation. The cylinder has a radius of 10 body lengths, which is the location of the most distant wind 
tunnel data that is available. The cylinder is clipped ahead of 32.6 body lengths aft of the model and behind 
41.0 body lengths aft of the model. The cylinder is restricted to its lower quadrant. The initial integration 
surface is poorly resolved due to the initial coarse grid. The pressure signature is not visibly propagated to 
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this initial integration surface. The final adapted grid integration surface is a much better approximation of 
a cylinder due to the background grid refinement. The peak signature pressure is larger at the horizon than 
the centerline because the model is designed to have a reduced centerline pressure signature. 




(a) Initial grid. (b) Final grid. 

Figure 20. SLSLE pressure on quarter cylinder integration surface 10 body lengths below the model. 


The adaptation history of the pressure integral and its error estimation is shown in Fig. 21. The requested 
error tolerance tolo is set to half of In at each adaptive iteration for grids sized less than 2,000,000. Above 
2,000,000 control volumes, the tola is set to /a, reducing the rate of grid growth. The goal of increasing tola 
is to obtain more resolved results with a more efficiently distributed and aligned grid at the expense of more 
adaptation cycles. 31 The change in requested error tolerance is observed as a reduction in grid growth per 
adaptive iteration in Fig. 21(a). This case shows a less dramatic reduction in the remaining error estimate 
over the final few grids than the previous cases, which may be due to the less aggressive tola = In on the 
final grids. 
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(a) Pressure integral. 


(b) Pressure integral remaining error. 


Figure 21. SLSLE pressure integral and uncertainty convergence at 10 body lengths. 


The pressure signatures at one body length for the series of grids employed during adaptation is shown 
in Fig. 22. One body length is much closer than the integration surface, but the signal must be resolved at 
this location for it to propagate to the integration surface. The pressure signature is broadly smeared on the 
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initial grid. The signals grow in amplitude with adaptation and the pressure peaks sharpen on the final few 
grids. The heuristic limiter is employed during adaptation. The signature mismatch near x/l = 0.8 is also 
present at other propagation distances so will be discussed with Fig. 24. 



Figure 22. SLSLE pressure signature adaptation history at 1.0 body lengths. 


The Venkatakrishnan, heuristic, and Barth-Jespersen limiters are applied on the final grid. The resulting 
pressure signatures at one body length are shown in Fig. 23. The signatures of all three limiters are very 
similar, except near discontinuities. The Barth-Jespersen limiter reduces the over- and under-shoots of the 
bow and tail shocks. The difference between the different limiter signatures is greater for this case than the 
cone cylinder and delta wing body cases. 



Figure 23. SLSLE final adapted pressure signature at 1.0 body lengths for various limiters. 


Centerline pressure signatures are presented in Fig. 24 for 1.0, 1.5, 2.0, and 2.5 body lengths below the 
model. Langley wind tunnel data is available for all four locations, but the closest Glenn wind tunnel data 
is available at 2.5 body lengths below the model. The Langley and Glenn wind tunnel measurements are 
generally in good agreement at 2.5 body lengths in Fig. 24, but the small differences of the two measurements 
gives an indication of the level of uncertainty in the measurements. The agreement between the wind tunnel 
and computed signatures is good, except in the region near x/l = 0.8. Other investigators 11,34 also showed 
a difference between wind tunnel and computed pressure signatures at x/l = 0.8 as well as provided posible 
explanations. The model geometry in this mismatch region is examined in Park, 37 which illustrated a 
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sensitivity in the signature near an uncertainty in the “as designed” CFD geometry and the “as built' 
geometry of the wind tunnel model. 




Figure 24. SLSLE centerline pressure signatures for various locations below the model. 


Figure 25 compares the adapted cut-cell method with the Glenn wind tunnel measurement at 10 body 
lengths. This is the same distance as the integration surface. The shock strength increases away from the 
configuration centerline. The front portion of the signature is well predicted in both Fig. 24 and 25. The 
aft portion of the computed signal shows the largest difference from the wind tunnel data at all propagation 
distances. The discrepancy between the wind tunnel and computed signatures near x/l = 0.8 decreases for 
the signatures away from the model centerline. The “as built” and “as designed” geometries correspond 
better laterally. 


VII. Summary and Concluding Remarks 

A cut-cell approach to CFD that utilizes the median dual of tetrahedral background grid with a topolog- 
ically consistent cut-cell geometry calculation method is described. The modifications made to the FUN3D 
inviscicl residual to account for the irregular geometry of the cut cells is detailed. The implications of re- 
construction limiting on the discrete adjoint solution are discussed. A single grid output adaptation scheme 
is described that eliminates the need to generate a memory intensive embedded grid. Output-based adap- 
tation simulations are presented for a cone-cylinder, delta wing, and low-boom configurations to predict 
off-body pressure signatures in supersonic inviscid flow. These predicted signatures are compared to wind 
tunnel measurements up to 10 body lengths below the model, providing simultaneous on- and off-centerline 
predictions. Multiple limiting schemes are applied to illustrate their effects. 

The use of cut-cells with an output-based adaptive scheme completely automates this accurate predic- 
tion capability after the triangular surface mesh is generated. The use of a heuristic limiter and cut cells 
with an anisotropically adapted background grid has dramatically improved the robustness of output-based 
adaptation as applied to sonic boom prediction. This allows extremely coarse isotropic initial grids that can 
be generated without a priori knowledge of shock locations or Mach angles. The general anisotropy of the 
adapted background grids allows for centerline and off-centerline signal prediction. This is the first published 
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Adapted Cut Cell - 

Glenn Wind Tunnel, h/l=10, 0.0 deg Off Center Line 




Glenn Wind Tunnel, h/l=10, 17.6 deg Off Center Line 



Glenn Wind Tunnel, h/l=10, 31.7 deg C 




-0.5 0 0.5 1 


Figure 25. SLSLE pressure signatures for various locations 10 body lengths below the model. 
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CFD prediction of the off-centerline SLSLE signatures at 10 body lengths. 
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